NASA Contractor Report 198392 


Modeling of Bulk Evaporation 
and Condensation 


S. Anghaie and Z. Ding 
University of Florida 
Gainesville, Florida 


April 1996 


Prepared for 

Lewis Research Center 

Under Contract NAS 3-263 14 



National Aeronautics and 
Space Administration 


TABLE OF CONTENTS 


Abstract 2 

Nomenclature 3 

CHAPTERS 

I. INTRODUCTION 5 

II. FORMULATION 9 

III. ANALYSIS AND ASSUMPTIONS 13 

IV. VAPOR PHASE FRACTION UPDATE 18 

V. RESULTS AND DISCUSSIONS 22 

5.1 Accuracy Issues 22 

5.2 Comparison of H-Formulation and E-Formulation 24 

5.3 Comparison of T-Based and E-Based Update Methods 24 

5.4 Bulk Evaporation, Condensation and the Interface Motion Scale 

Analysis 28 

VI. CONCLUSIONS 31 

Acknowledgment 33 

REFERENCES 34 


1 



ABSTRACT 


This report describes the modeling and mathematical formulation of the bulk evaporation 
and condensation involved in liquid-vapor phase change processes. An internal energy 
formulation, for these phase change processes that occur under the constraint of constant volume, 
was studied. Compared to the enthalpy formulation, the internal energy formulation has a more 
concise and compact form. The velocity and time scales of the interface movement were obtained 
through scaling analysis and verified by performing detailed numerical experiments. The 
convection effect induced by the density change was analyzed and found to be negligible 
compared to the conduction effect. Two iterative methods for updating the value of the vapor 
phase fraction, the energy based (E-based) and temperature based (T-based) methods, were 
investigated. Numerical experiments revealed that for the evaporation and condensation problems 
the E-based method is superior to the T-based method in terms of computational efficiency The 
internal energy formulation and the E-based method were used to compute the bulk evaporation 
and condensation processes under different conditions. The evolution of the phase change 
processes was investigated. This work provided a basis for the modeling of thermal performance 
of multi-phase nuclear fuel elements under variable gravity conditions, in which the buoyancy 
convection due to gravity effects and internal heating are involved. 
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NOMENCLATURE 


A, C,D 

E 

H 

O 

T 

V 

U 

W 

X, R 



h 


f 

l 

P 

t 

x, r 


constants 
internal energy [J] 
enthalpy [J] 
heat quantity [J] 
temperature [K] 
volume [m 3 ] 
velocity [m/s] 
work [J] 

dimensionless coordinates 

specific heat at constant volume and constant pressure [J/kg/K] 

specific internal energy [J/kg] 

specific enthalpy [J/kg] 

vapor phase fraction 

length scale [m] 

pressure [bar] 

time [s] 

cylindrical coordinates [m] 


Greek symbols 

a thermal diffusivity [m 2 /s], or thermal diffusivity ratio 


phase-change temperature interval 
thermal conductivity [W/K/m], or its ratio 
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0 
0 

p 

x 

Subscripts 

/ 

1 

l 

r 

s, sat 
v, vap 

Superscripts 

k 

n 


order of magnitude 
dimensionless temperature 
density [kg/m 3 ], or its ratio 
dimensionless time 

interface position, mixed phase 

interface 

liquid phase 

characteristic scale 

saturation point 

vapor phase 

iteration level 
time level 
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CHAPTER 1 


INTRODUCTION 


In recent years, phase change processes such as melting and solidification, boiling and 
condensation have attracted increasing interest in their development and applications Among 
these activities, numerical modeling and simulation of the phase change processes have been the 
focus of many investigators. The purpose of this work is to develop a consistent numerical 
method for the computation of the bulk evaporation and condensation processes involved in a 
specific class of liquid-vapor phase change problems under constant volume and no buoyancy 
conditions. The need for the development of such models has been promoted by a recent 
investigation aiming to establish the feasibility of using an encapsulated two-phase nuclear fuel for 
high temperature space power and propulsion applications. 

Numerical modeling of heat transfer problems involving phase. change has been discussed 
by a large number of investigators [1-14]. Shyy [11] presented a comprehensive discussion of the 
various computational elements important for the prediction of complex fluid flows and interfacial 
transport. Most of the noted investigators have approached the solid-liquid phase change 
problems, such as those involved in melting and solidification processes. There has also been an 
ongoing interest in the investigation of liquid-vapor phase change problems, such as evaporation 
and condensation More attempts have been made to numerically model and simulate 
condensation phenomena In particular, a number of investigators have developed numerical 
models for flow condensation problems [15-18], Most of this numerical modeling has focused on 
parabolic flow condensation problems where the forced convection plays a dominant role. In this 
type of condensation problem, the vapor pressure, saturation temperature and latent heat have 
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been treated as constant properties. In the present study, the saturation temperature, and pressure 
are modeled as dependable variables. In a more recent work, Jang et al. [19] have modeled a heat 
pipe start-up from the frozen state by using the a system of simplified one dimensional equations. 
Shyy et al. [20] have successfully developed a computational method for predicting the two-phase 
transient flow and heat transfer characteristics within a constant pressure reservoir of a capillary- 
pumped-loop and investigated the phase change of liquid and vapor, under various conditions. 
The present approach has expanded the mathematical modeling of the phase change problems by 
including the pressure as a dependent variable and using the internal energy formulation. 

The numerical models for the phase change problems can be divided into two groups, 
based on the choice of dependent variables. In the first group, temperature is the sole dependent 
variable, and energy equations are written separately for the solid and liquid phases. In the second 
group, enthalpy is used as the main dependent variable along with the temperature. This 
formulation which is widely used by many investigators is known as the enthalpy method. In the 
enthalpy model, the phasic interface is eliminated from the formulation and the problem is made 
equivalent to a nonlinear heat conduction problem without phase change. Shamsunder and 
Sparrow [3] and Shyy [11] summarized these features and demonstrated the usefulness of the 
enthalpy formulation. The enthalpy formulation has been employed in many investigations [3, 9, 
10 , 11 , 20 ], 

This report presents a discussion of the bulk evaporation and condensation processes 
taking place in a constant volume container without vapor bubble or liquid droplet generation and 
boiling. This problem can be depicted with the case presented in Fig. 1 . A rigid, cylindrical 
container is partially filled with vapor and liquid. Initially, the vapor and liquid are kept saturated 
at a certain pressure and uniform temperature. Well selected boundary conditions on the walls are 
imposed to activate the pure bulk evaporation or condensation, and also to avoid the buoyancy 
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Figure 1 . 1 Schematic of the bulk evaporation and condensation 

effect and bubble generation. While keeping the side wall insulated, only heating the top wall or 
cooling the bottom wall induces pure evaporation or condensation, respectively. The volumetric 
expansion or contraction due to the density change between two phases induces convection in the 
vapor phase. The scaling analysis conducted in this work indicates that this convection effect is 
quite week and negligible. It can be reasonably assumed that the phase change processes are 
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driven entirely by the conduction and the release or absorption of latent heat at the phase change 
interface. Similar situations are also possible for the multi-dimentional phase-change problems in 
zero gravity condition. 

No expansion work is involved in this constant volume phase change process. In the 
absence of mechanical work, the first law of thermodynamics indicates that the heat exchange 
between the system and the ambient equals the change of the internal energy of the system. The 
released or absorbed latent heat associated with the phase-change is equal to the change of 
internal energy. Based on the physical conditions of this system, a compact internal energy 
formulation of the energy equation is proposed. This formulation, similar to the enthalpy 
formulation, is a single region formulation, where one set of governing equations can describe 
both phases. 
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CHAPTER 2 


FORMULATION 


The first law of thermodynamics, for any thermodynamic system, is 

SQ = dE+5W = dH-(pdV + Vdp)+SW (2.1) 

where E is the internal energy, H is enthalpy, 0 is the heat input, W is the work output, p and V 
are the pressure and volume of the system, respectively. For the phase change process involved in 
the rigid and enclosed container with constant volume, there is no mechanical work. Therefore the 
first law for this closed and rigid system reduces to 

SO = dE = dH - V dp (2.2) 

which indicates that the heat exchange between the system and the ambient causes the internal 
energy, the enthalpy and pressure of the system to change. Therefore the released or absorbed 
latent heat associated with the phase change induces the changes of internal energy, enthalpy and 
pressure. The quantity of the heat exchange is equal to the internal energy change, but not equal 
to the enthalpy change. The conventional enthalpy formulation of the energy equation, commonly 
used for the phase change between solid and liquid, should include a transient pressure term for 
this phase change between liquid and vapor. The weak form of the energy equation based on 
enthalpy formulation is 

Dh = Dp + V n (2.3) 

y Dt Dt 

This formulation including the transient pressure term introduces additional complexity 
and does not seem to be the best choice for this problem. A more attractive alternative, the 
internal energy based formulation for the energy equation which does not include the pressure 
term, is proposed and developed in this work. Similar to the enthalpy formulation used for melting 
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and solidification problems, the internal energy formulation does not require to explicitly track the 
interface position. 

The weak form of the energy equation with internal energy as the primary variable can be 
expressed as follows: 

P—=V(kVT) (2.4) 

Dt 

where the internal energy e may be expressed, in the form of a summation of sensible heat and 
latent heat, for the liquid, vapor and mixed phases as: 

e = c vl T / = 0 

e = c vJ T sat + f Ae 0 </< ( 2 . 5 ) 

e = c T + Ae + c (T-T ) / = 1 

v, / sat v , v sat 

where T sat is the saturation temperature, Ae is the internal energy difference between saturated 
vapor and saturated liquid (latent heat), and the constant volume specific heat c v ’s are functions of 
temperature provided by Reynolds [21], The vapor phase fraction /is zero in the region filled with 
liquid phase and unity in the region occupied by vapor phase. The vapor phase fraction lies 
between zero and unity when the control volume is undergoing phase change. 

Substituting Eq.(1.5) into the governing Eq.(1.4) yields 


= ?■(*, vr> 

f = o 


P,^.r T ) = V(K f VD-p f lxSL 

0 < / < 1 

(2.6) 

P,j^(c r ..T) = v(K.vr) 

/= 1 



where the latent heat now appears as a source term. The form of this internal energy based 
equation is identical to that of the enthalpy based equation, except the pressure term appearing in 
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the later one. Just like the enthalpy formulation, this internal energy formulation is also a single 
region formulation, where one set of governing equations can describe both phases. A major 
advantage of the internal energy formulation is that it is not necessary to split the domain into 
separate subdomains consisting of different phases. A fixed grid can be employed to facilitate the 
computations. 

Both enthalpy formulation and internal energy formulation are implemented and examined 
The vapor phase fraction, f may be defined in two different ways, which can lead to different 
update methods. 

Some additional closure relations are required for the multiphase condition in the 
container. Property variations are based on a simplified searching routine. The different phases are 
differentiated based on the local temperature or internal energy at each nodal point, within the 
frame work of the temperature based or internal energy based update method. Using temperature 


based update, the liquid or the vapor phase is identified by the local temperature, wherever the 
temperature at any point is lower than the saturated liquid temperature or higher than the 
saturated vapor temperature, that point is located in the liquid or vapor region, respectively. 
Similarly, if the internal energy based update is used, the liquid or the vapor phase is identified by 
the local internal energy. Either case allows the vapor phase fraction in a computational cell to be 
calculated. 

Using the Clausius-Clapeyron equation, the saturation temperature as a function of the 
system pressure is calculated. 


C 

Tsat = — ' 
D-\np sat 


( 2 . 7 ) 


where C and D are constants. 
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The saturation pressure, or the vapor pressure, can be obtained from the approximate 
equation of state for an ideal gas, and the mass conservation in which the bulk vapor density and 
temperature are used: 

Psat = P™p = m tota i = m v + m/ = Constant (2.8) 

Y vap 

where R is the gas constant for vapor phase. Eqs.(2.7) and (2.8), have been proposed and used as 
the basis for the model developed by the authors, Ding and Anghaie [ 22, 23 ]. 
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CHAPTER 3 ANALYSIS AND DISCUSSION 

To simplify and generalize the formulation for the bulk evaporation and condensation 
processes, these governing equations and closure relationships should be normalized Scaling and 
analysis of order of magnitude are conducted to develop simplifying assumptions that lead to a 
concise and efficient solution procedure for the governing equations. For the different phases of 
any material, the thermodynamic and transport properties are quite different, and they may be in 
the same order or may have a difference of several orders of magnitude. For the same phase, some 
properties are independent or slightly dependent on temperature or pressure, but others may be 
strongly dependent. The following relations are found about the orders and variations of these 
properties, for liquid and vapor phases far from the thermodynamic critical point, within a 
moderate phase-change range. 


c v , = constant , 

c v v = cons tan?. 

C V / 

— 0(1), 

c v,v 

K t = constant , 

k v = cons tan / , 

*7 

K = 

~0(10)-0(10 2 ), 

p; = cons tan f. 

p v * constant. 

II 

< i*-- 

~0(10 2 )-0(10 3 ), 

a , = constant , 

a v * constant. 

a = 

~ O(10 _1 ) 


a v 


where ‘ = constant ' means that the variation is within one order of magnitude, and ‘ O 
represents the order of magnitude. 

Therefore all the properties, except the vapor density and vapor thermal diffusivity, can be 
assumed constant. The liquid specific heat and vapor specific heat are assumed equal. In the 

following analysis and computation, these values of the ratios such as 

K _ El = \Q,.p 0 = — = 10 2 ,.a 0 = = 10~‘ are used, where the subscript 0 denotes a 

K v P v fl ^v.O 

reference state 

As discussed in some literature [1, 24], due to the density difference between two 
different phases, the volumetric expansion effect induces the convection in the phase of lower 
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density. Let Ui denote the speed for the motion of the interface and Uv for the motion of the 
vapor phase resulting from the density change. The mass conservation at the interface should be 
satisfied as follows 

A U,=PM-U.) (3.2) 


The energy-balance equation at the interface is 


OF, OT v T J ( \TT 

K l~T ~ K ^ 7 + Pv e v U v = (P*«* “ Pi e i) U , 
Ox Ox 


( 3 . 3 ) 


where e v and denote the internal energies for the saturated vapor and liquid phases, 
respectively. 

Eliminating U v from Eq.(3.3) via Eq.(3.2) yields 


fc.v QlJ ~ 


= K 


oT v a j t 

= PA eU i 

ox ox 


( 3 . 4 ) 


where q,, v and q,j represents the heat flux in the interface at the vapor side and liquid side, 


respectively. 

Eq.(3.4) can be rearranged as 


Pl&e ox ox l t aX l t p OK 


( 3 . 5 ) 


where /, is the thermal diffusion length scale which is the length of the container, St=c v T</As, 


9=T/To, X=x/l t . p=p/p v . There are two possible velocity scales at the interface, i.e., Ui=St(a t! /lJ, 


U 2 =St(aVh)/p , and their ratio is — = = a p ~ 0(10 ')-O(10 2 ) ~ 0(10) , Ui>U 2 

U 2 a v p v 

Therefore it is reasonable to assume that the lower velocity scale U 2 controls the movement of 
the interface. 
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Eq.(3.5) can also be expressed as 


u =U ] ^ l -U 2 ^ l 
' 1 dX 2 dX 


(3.6) 


The velocity, U 2 =St(aJl<)/p, includes the density change effect which actually reveals the 
constraints of energy balance and mass conservation at the interface, and the volumetric 
expansion or contraction during evaporation or condensation, respectively. It is lower than the 
pure thermal diffusion velocity in liquid, a/l t and that in vapor a// f . Thus U 2 is to be chosen as 
the characteristic velocity scale where the density change effect, as well as the S, number, appears 
explicitly. The dimensionless velocity of the interface movement is 




1 U 2 d X dX 


(3.7) 


For the case of evaporation, for example, heat is supplied to the system from the top wall 
which is the boundary for the vapor phase. If the appropriate amount of heat flux is provided to 

d6„ 


ensure 


ax 


0(1), it is found that U, ~ 0(1), thus Uz ~ Ur = S t (cc v / /, ) / p Numerical results 


which will be shown in the next chapter are compatible with this expression. 
The characteristic time scale is 

t 

' U r cc v S t 


(3.8) 


Obviously this time scale is larger than those of thermal diffusion in vapor and liquid 


I* j 2 j 2 

phases that is — p> — > — , which implies that the motion of the interface controls the rate 
H a.. 




of tjijs heat transfer process. From the above analysis, it can be conclude that the motion of the 
interface is affected by the density difference between liquid and vapor phases, thus its velocity is 
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very low. During the evaporation process with constant system volume, density change increases 
the vapor pressure and the saturation temperature, which in turn suppresses the evaporation. 

From Eqs.(3.2) and (3.7), the velocity of vapor convection induced by the density change 

is Uv - (1 - p)Ui ~ - — ——S t , and since p »1, U v ~ (a^/lJSt , which is still very low. 

P h 

Therefore the Reynolds number and Peclet number for the convection are Re= U v l v /v v =(iyi)(St 
IPr), Pe= (l/l) St. , where l v is the vapor motion length scale. Since the velocity of the vapor 
convection is very low, and the convection just occurs near the interface, the vapor convection 
length scale will be very small compared to the thermal diffusion length scale, say, ///, < 0(1). 
For most vapors, Pr ~ 0(1), and if moderate St is considered, i.e., St ~ 0(1), then Re < 0(1), 
and Pe < 0(1). For the problems in which Re and Pe are less than the order of unity, the 
convection effect is weaker than the molecular movement effect and hence insignificant. Therefore 
the thermal diffusion, i.e., the heat conduction, dominates the heat transfer process and the 
convection effect induced by the density change is negligible. 

On the other hand, it can be found that the convection effect caused by the density change 
during phase change is much less significant than the convection induced by buoyancy effect when 
encountered. The characteristic velocity of the natural convection is Ur=(Gr) 1/2 v/L where Gr is 
the Grashof number, the corresponding Reynolds number Re=UrL/v =(Gr) 1/2 . If Gr=10 4 ~10 6 , 
then Re^ltf-lO 3 , which is stronger than the density change effect. 

From this analysis, the convection effect can be reasonably neglected and the bulk 
evaporation and condensation processes can be assumed entirely driven by heat conduction. The 
governing equations (2.6) can be reduced to the heat conduction-driven phase-change problems: 
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Eq.(3.9) can also be written in terms of dimensionless parameters as . 


S t dd 2 
— =aV L e 

p &z 

EsJL.'Ja-L* 1 

p 8x p fa 

p ck 


f = o 

0</<l 
/= 1 


(3.10) 


Here S, is the Stefan number defined as c v T s yAs , where T s ,o is the initial saturation 
temperature. T s , 0 and the length of the container are chosen as the characteristic temperature 
scale and length scale, respectively. The characteristic time scale is shown in Eq.(16). The 
relationships presented in Eqs. (2.7) and (2.8) are nondimensionalized as 




A-\n(p s ) 


Ps = Pv e 


(3.11) 


where p s and p v are dimensionless saturation pressure and vapor density, respectively. 



CHAPTER 4 


VAPOR PHASE FRACTION UPDATE 


To solve the physical model described by Eq.(3.10), a closure relationship between the 
vapor phase fraction / and the temperature 9 or the internal energy e should be formulated. For a 
pure substance undergoing evaporation or condensation, the total internal energy is a 
discontinuous function of the temperature: 


e = — = S t - 9 

Ae 

e — — — S t • 9 + \ 
Ae 


6 < 9 sat 


(4.1) 


e>e 


sat 


where 9 sat is the saturation temperature. However, from a computational viewpoint, 
discontinuities are difficult to track and it is often necessary to smear the phase change over a 
small temperature range to attain numerical stability. Thus we can have 


e = S t -9 Q<Qi 

e = S t -9 l +(9-9 l )/(9 v -9 1 ) 0/<0<0 v (4.2) 

e = S t -9 + 1 6 V <8 

where 0/= 9 sa r£, 9 v =9 sat +£, and £ is the phase change interval. As a result, the discontinuity is 
replaced by a small window, 2s, over which the phase change occurs. For a pure substance, this 
window is a purely numerical artifact and must be chosen as small as possible in order to 
accurately model the physical system. Substitution of Eq (4.2) into the normalized form of 
Eq.(2.5) yields 
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e-e l e-e l 


(4.3) 


/ = 

e v - e t 2 s 

which is used to iteratively update the phase fraction from the computed temperature field. This 
technique has been successfully used by Shyy et al. [10, 11, 20] and was referred to as the T- 
based update method This technique has proven to be very effective for achieving convergence 
An alternative formulation to the T-based update method is to express the phase fraction as a 
function of the total enthalpy rather than the temperature. This update is referred to as the H- 
based update method. More information about this method can be found in Ref. [10, 11] 

A third formulation expressing the vapor phase fraction as a function of the total internal 
energy, rather than the temperature or enthalpy, is developed in this work. According to the 
classical thermodynamic definition, the relationship between the vapor phase fraction and the 
internal energy is 

f = e-e l (4.4) 

This is used for the iterative update of the vapor phase fraction This update procedure is 
referred to as the E-based update method and is similar to the isothermal case discussed by Voller 
and Swaminathan [2], In particular, it is noted that e=0 is allowable in this formulation, thereby 
allowing for modeling phase change of a pure substance with higher accuracy. This is the most 
important advantage of the E-based method compared to the T-based method where a significant 
effort must be made to determine the phase change temperature interval s. For the phase change 
problems involving complex interface geometry, the normal temperature gradient varies with the 
positions along the interface The determination of the phase change intervals for such problems is 
a tedious task. The E-based method does not need to use a phase change window and can be 
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conveniently used for all kinds of phase change problems including the ones with complicated 
geometry. Another advantage of the E-based method is the fact that if the latent heat is large 
relative to the sensible internal energy, the vapor phase fraction varies slowly with respect to the 
total internal energy. This enhances the stability of the update procedure. For the liquid-vapor 
phase change occurring far away from its critical point, the latent heat is relatively large and the 
application of the E-based update method even for simple geometry is recommended. In the 
present work, the T-based and E-based methods are examined in details from the viewpoint of 
accuracy and computational efficiency. 

A straight-forward implementation of the T-based method to discretize Eq. (3.10) for 
each computational cell yields: 

+i* (4 5) 

p Ar 

where the subscript P represents the nodal point of the control volume to be solved, and nb 
represents the neiboring nodal points. This equation is solved according to the T-based update 
procedure given in Eq.(4.3). As reported by Shyy et al.[9, 10] and Prakash et al. [5], for Stefan 
number less than unity (large values of latent heat), this form of the iterative update is prone to 
numerical instability in the form of oscillations. For moderate values of Stefan numbers, the use of 
an under-relaxation approach for the calculation of the source term facilitates the solution 
process. The main cause of the oscillatory, non-convergent behavior is the extreme sensitivity of 
the vapor fraction to small changes in temperature which leads to a large negative feedback effect 
through the source term. 
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Substituting Eq.(22) into the discrete form of the conservation law and moving the 
resulting coefficients of T p into the a P term in Eq.(24), yields a set of stable equations similar to 
those formulated by Shyy et al. [9], These equations are: 


a P 6 p k = 1**0 J +b k - 


\cip + 


A V r 

r p 


P&T 


AV 1 
pAr 2s 



< 6\,or6p 

h ^r [ 7 L+// ’ ,] 

pAr 2s 

VI 

7 

VI 


(4.6) 


where the superscript n indicates the current time step and k indicates the current iterative value 
Using this form of the update procedure, the temperature and the vapor phase fraction fields can 
be computed simultaneously and are therefore coupled at every step of the iterative process. This 
numerical procedure converges quickly with no loss of accuracy, granted that a proper value of s 


is used. 


The implementation of the E-based method is relatively straightforward and does not need 
any pre-conditioning. The internal energy is updated first, with the current values of saturation 
temperature and latent heat, 

e p = S,6 P + /*-' , 0 < f k ~ l < (4.7) 

Then the vapor phase fraction is updated with the current internal energy, 

f k = e p - e k , e k < e p < e k (4.8) 

where ei=St6 sa , and e^epAe. 
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CHAPTER 5 


RESULTS AND DISCUSSION 


5.1 Accuracy Issues 

A testing case of pure bulk evaporation is selected to examine the accuracy, the 
computational performances of each formulation with different update methods. A fixed volume 
container with one fourth volume filled with a saturated liquid is analyzed. The rest of the 
container’s volume is filled with the saturated vapor. The initial temperature of the saturated fluid 
is 1.0. The top wall is heated and kept at a temperature of 1.5, while the bottom wall is insulated. 
Such an orientation of heat flow is required to avoid the bubble generation and boiling, while the 
bulk evaporation is allowed. A mesh of 81x6 non-uniform and fixed grids is adopted for the liquid 
and vapor within the container. The time step size Ax=10' 4 is chosen. The second order central 
differences are used for all spatial derivative terms and the fully implicit difference is employed for 
the time marching. 

Figs. 5. 1(a) and (b) show the evolution of the evaporation process. At the time instant of 
0.03, about 18% of liquid is vaporized, and the saturation temperature is about 1 .38. A little wavy 
or non-uniform variation can be seen due to the nature of the discretization and the non-uniform 
release or absorption of the latent heat. A significant part of the latent heat content in a 
computational cell is absorbed when it starts to undergo evaporation. Thus the interface 
movement is slowed down when phase change in this cell starts and subsequently speeds up 

The grids are clustered near the liquid-vapor interface, and the finer grids are placed near 
the interface and covering the range of the interface movement. Different minimum grid sizes, i.e., 
AX m ,„= 1/100, 1/200, 1/400, are selected to examine the solution accuracy and the independence 
of solution on the grid size. The results with different minimum grid sizes are basically identical, 
while the finer grid size yields smoother solution. In the following computations, the grid size with 
AX mm =1/200 is used. 
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Saturation Temperature Evaporated Liquid Fraction 



Time (x) 



Time (x) 


Figure 5. 1 Solution features for bulk evaporation (Grids 81x8, AX min -l/120, 1/200, 1/400) 




5.2 Comparison of H-Formulation and E-Formulation 


The same testing case is computed to examine the performances of different formulations, 
the enthalpy formulation and internal energy formulation. Figs. 5. 2(a) and (b) show the evaporated 
liquid fraction and the variation of the saturation temperature during the evaporation process. 
Fig.5.2(c) shows the efficiency of the computations, that is, the number of iterations required at 
each time step for the residual to go below lCf 4 . The enthalpy formulation needs a few more 
iterations for convergence of the calculation. It can be found that the internal energy formulation 
and enthalpy formulation yield nearly identical results, with close computational efficiency. It 
indicates that both of these formulations are proper for this problem and have similar efficiency 
performance. The only advantage in the internal energy formulation is the more concise form 
without an explicit pressure term In the following investigations, only the internal energy 
formulation is adopted. 


5.3 Comparison of E-Based and T-Based Update Methods 

The T-based update method and the E-based update method are examined for the same 
testing case. Figures 5.3(a) and (b) show that the T-based method and E-based method yield 
nearly identical results. It can be seen that for convergence the E-based method requires about 
half of the iterations which is needed for the T-based method. Thus the E-based method is 
computationally more efficient than the T-based method, as shown in Fig. 5. 3(c). The efficiency 
performances of these two update methods for the evaporation and condensation problems are 
different from those for melting and solidification problems where the T-based is better than H- 
based method as reported [10]. 
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The reason for the improved efficiency of the E-based method is that the phase change 
temperature window is used for the T-based method, which is not required for the E-based 
method. When the T-based method is used and the purely artificial phase change window has to 
be adopted, it is not easy to use the temperature to identify the phase interface since the 
temperature is always continuous at the interface For melting and solidification problems, the 
phase change temperature is usually constant. The saturation temperature, however, changes 
continuously for evaporation and condensation processes in a constant volume system. This 
difference causes distinctive computational efficiency performances. For each iteration at every 
time step, the saturated vapor and liquid temperature as well as the saturation temperature should 
be updated Therefore, the T-based method converges more slowly. When the E-based method is 
used, the internal energy has a sharp jump at the interface, it is easy to use internal energy to 
identify the phases As has been mentioned before, the relative large amount of the latent heat will 
result in a more stable update procedure and faster convergence rate. Furthermore, as discussed in 
the section of update methods, if the T-based method is used, there will be some difficulty to 
determine the phase change temperature intervals for the complicated phase change problems 
where the interface could be of any arbitrary shape. Using the E-based method- can completely 
circumvent the need for a phase change window. Therefore, the E-based method is superior to the 
T-based method, in terms of computational efficiency and formulation simplicity, and is capable of 
handling more complicated phase change geometry 
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(a) evolution of vaporization 



(b) variation of saturation temperature 



(c) computation efficiency performance 

Figure 5.2 Comparison of two formulations for evaporation 
The solid line: the internal energy formulation, the dashed line: the enthalpy formulation 
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(a) evolution of vaporization 



(b) variation of saturation temperature 



(c) computation efficiency performance 

Figure 5.3 Comparison of two update methods for evaporation 
The solid line: the E-based method, the dashed line the T-based method 
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5.4 Bulk Evaporation, Condensation and the Interface Motion Scale Analysis 

In the following work, the internal energy formulation and the E-based update method are 
used to investigate the bulk evaporation and condensation processes and also to analyze the 
interface motion scale. With other conditions remaining the same, i.e., the initial temperature, the 
grid distribution and the time step size, different boundary conditions are imposed at the top and 

bottom walls. First, the container’s top wall is heated with a constant heat flux 

dO 

a = — -I „ = 1 0 and the bottom wall is insulated. After the time instant of 0.02, the top wall 

~ tw \ <*=U 7 

is insulated while the bottom wall is cooled with a constant heat flux 

q bw =,c— L 10 = 10x1.0= 10. Correspondingly, the system undergoes an evaporation- 
dX 

condensation process. Fig. 5. 4(a) shows the variation of the saturation temperature during the 
transient. The motion of the interface, which is downward and upward during the evaporation and 
condensation, respectively, is exhibited in Fig.5.4(b). The axial temperature distributions at 
different time instants during the evaporation process are plotted in Fig. 5. 4(c), and those of 
condensation process are plotted in Fig.5.4(d). 

The velocity of the interface motion can be calculated from results presented in Fig . 5.4(b). 
During the evaporation and condensation processes, the interface velocity is always constant 
because only the constant heat flux at either the top wall or bottom wall is provided. During the 
evaporation process which takes place with a constant heat flux at the top wall, the interface 
velocity is a constant about 2.5, which is in the order of 1.0. While during the condensation 
process with a constant heat flux at the bottom wall, the interface velocity still remains constant at 
the level of about 0.8. This indicates clearly that the nondimensional interface velocity is 
consistently of the order of 1.0. In other words, it can be concluded 
that U, ~ 0(1), thus U i ~ Ur = S,(a v / /,) / p. This analysis verifies the phenomenological 

analysis presented in previous chapter of this report. This also indicates that the velocity 
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scale Ur = S t (a v 1 1,)/ p has been chosen appropriately so that the velocity of the interface 
motion can be calculated directly. 
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(a) Variation of saturation temperature during evaporation and condensation 
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(b) Motion of interface during evaporation and condensation 



Figure 5.4 Bulk evaporation and condensation 
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Temperature (0) Temperature (0) 
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(c) Temperature distribution during evaporation 



Axial position X 


(d) Temperature distribution during condensation 


Figure 5.4 ( continued ) Bulk evaporation and condensation 




CHAPTER 6 


CONCLUSIONS 


In this work, a numerical method was developed to model the bulk evaporation and 
condensation processes in an encapsulated container. An internal energy formulation for the 
constant volume phase change process was successfully implemented. Some assumptions 
regarding the convection heat transfer rate at the vapor-liquid interface were made through 
scaling and order analysis The velocity and time scales of the interface movement were obtained 
through scaling analysis. The dominating mechanism for the heat transfer at the interface is 
conduction which is modeled in the dimensionless governing equations. The convection effect 
induced by the density change was analyzed and found negligible. 

The computational performances of the well-established enthalpy formulation and the 
proposed internal energy formulation were evaluated and compared. Both formulations yield 
identical results with similar computational efficiencies. The internal energy formulation has the 
advantage of having more compact form. With the internal energy formulation, the internal energy 
based update method was developed The performances of two update methods for the vapor 
phase fraction, the E-based and T-based methods, were investigated. The E-based method is 
superior to the T-based method in terms of computational efficiency and formulation simplicity 
for evaporation and condensation problems. The E-based method is also capable of handling more 
complicated phase change problems, because it does not require the artificial phase change 

window. 

The internal energy formulation and the E-based method were successfully implemented 
and used to compute bulk evaporation and condensation processes under different conditions. The 
evolution of the bulk evaporation and condensation processes was simulated. The variations of 
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the temperature distribution, the liquid content and the saturation temperature with time for a 
variety of boundary conditions were obtained. Numerical experiments and analysis verified the 
velocity scale of the interface motion. 

This work provided a basis for further modeling of other bulk liquid-vapor phase change 
problems of multi-phase nuclear fuel elements, in which the buoyancy convection due to gravity 
effects and internal heating are involved. The modeling of thermal performace of multi-phase 
nuclear fuel elements under variable gravity conditions is discussed in another report. 
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